#####################################################################
###### REPLICATION CODE FOR "CONVERGENT FLEXIBILITY: HOW ############
###### INTERNATIONAL LAW KEEPS PACE WITH TECHNOLOGICAL CHANGE #######
#####################################################################

## This file can be used to reproduce all figures/tables in the manuscript/appendix.

rm(list=ls())
graphics.off()
library(plyr)
library(dplyr)
library(tidyr)
library(broom)
library(ggplot2)
set.seed(888)  
options("scipen"=100, "digits"=4) 



######## PREPARING DATA #########

setwd("/Users/user/file location/") # change as needed


## Load file from registration survey (pretreatment)
pre <- read.csv("cf_registration.csv") 

## Load file for experiment 
post <- read.csv("cf_experiment.csv") 
intersect(names(pre), names(post))


## Covariates preceded by "X_"; separate and alphabetize
x_cols <- grep("^X_", names(pre), value = TRUE)
non_x_cols <- setdiff(names(pre), x_cols)
x_sorted <- sort(x_cols)
new_order <- c(non_x_cols, x_sorted)
pre <- pre[, new_order, drop = FALSE]
pre$X_LocationISO <- NULL


## Treat "" as NA across both files
pre <- pre %>%
  mutate(across(where(is.factor), as.character)) %>%                 # factors -> character
  mutate(across(where(is.character), ~ na_if(trimws(.x), ""))) %>%   # "" -> NA in character cols
  mutate(across(where(is.character), ~ na_if(.x, "NA")))             # optional: literal "NA" -> NA

post <- post %>%
  mutate(across(where(is.factor), as.character)) %>%                 # factors -> character
  mutate(across(where(is.character), ~ na_if(trimws(.x), ""))) %>%   # "" -> NA in character cols
  mutate(across(where(is.character), ~ na_if(.x, "NA")))             # optional: literal "NA" -> NA


## Merge into a single dataframe for easy analysis
df <- merge(pre, post, by=c("Synthetic_ID", "LawSocID"), all=T)


## Scale aversion scores (available for latter waves)
df$X_Aversion_Risk <- scale(df$X_Aversion_Risk, center = TRUE)
df$X_Aversion_Ambiguity <- scale(df$X_Aversion_Ambiguity, center = TRUE)

ncol(df)



####### PREVIEW IMPORTED DATA (pre-merged) ########

## Nyrup et al (BJPS, 2025): Cabinet officials since 1960 with terminal law degrees
hist(df$X_LawPtP)

## Polity V
hist(df$X_Polity)

## Guerriero (Data in Brief, 2016): Legal tradition
hist(df$X_ShariaLaw)
hist(df$X_CommonLaw)
hist(df$X_PureCivil)


####### ORGANIZATION OF REPLICATION FILE ########

# PAPER
# Table 2: Sample Statistics
# Figure 6: Point Range Plots
# Figure 7: See APPX 4.

# APPX 1
# A.1.2: Randomization Inference
# A.1.3: Balance Plots
# A.1.4: LASSO Regularization
# A.1.5: Numerical Results

# APPX 2
# Table A.2: Descriptive Statistics
# A.2.1: Wave Stability
# A.2.2: Nonparticipation 
# A.2.3: Manipulation Checks

# APPX 3 - Heterogeneous Effects

# APPX 4
# A.4.1: Topic Models
# A.4.2: Excerpts
# A.4.3: Mediation Analysis 




################################################
###### Tables 2 & A2: Treated Sample Stats ##### 
################################################

library(xtable)
library(purrr)
library(dplyr)

## Prepare data for table
dx_tbl <- df %>% dplyr::select(contains(c("X_", "Wave", "Treated")))

## Add helper variable
dx_tbl <- dx_tbl %>%
  mutate(
    Wave = case_when(
      grepl("^Pilot", Wave) ~ "Pilot",
      grepl("2020", Wave) ~ "Fall 2020",
      TRUE ~ "After 2020"
    )
  )

unique(dx_tbl$Wave)

## Prep wave order variable
dx_tbl$Wave <- gsub("Wave", "", dx_tbl$Wave)
dx_tbl <- fastDummies::dummy_cols(dx_tbl, select_columns = "Wave", remove_selected_columns = TRUE)

## Convert all columns to factors
dx_tbl[] <- lapply(dx_tbl, as.character)
dx_tbl[] <- lapply(dx_tbl, as.factor)

## Convert columns that do not contain any alphabetic characters to numeric
dx_tbl[] <- lapply(dx_tbl, function(x) {
  if(all(!grepl("[a-zA-Z]", levels(x)))) { # Check if any level contains alphabetic characters
    as.numeric(as.character(x))  # Convert to numeric if no alphabetic characters
  } else {
    x  # Keep as factor if there are alphabetic characters
  }
})

## Function to generate summaries for numeric and factor variables
summarize_data <- function(data) {
  data %>%
    dplyr::select(where(is.numeric)) %>%   # only numeric vars
    dplyr::summarise(across(
      .cols = everything(),
      .fns = list(
        N       = ~sum(!is.na(.)),
        Mean    = ~mean(., na.rm = TRUE),
        Median  = ~median(., na.rm = TRUE),
        `Std Dev` = ~sd(., na.rm = TRUE),
        Min     = ~min(., na.rm = TRUE),
        Max     = ~max(., na.rm = TRUE)
      ),
      .names = "{.col}_{.fn}"
    )) %>%
    tidyr::pivot_longer(
      cols = everything(),
      names_to = c("Variable", ".value"),
      names_pattern = "(.*)_(.*)"
    ) %>%
    dplyr::mutate(
      Variable = sub("^X_", "", Variable)
    )
}


## Apply the function on treated sample
dx_tbl_treated <- subset(dx_tbl, Treated==1)
dx_tbl_treated <- summarize_data(dx_tbl_treated)

rownames(dx_tbl_treated) <- NULL
xtable(dx_tbl_treated) %>%
  print(include.rownames = FALSE)

## As percentages
round(prop.table(table(subset(dx_tbl, Treated==1)$LawSocID)) * 100, 2)


## Plot table
print(xtable(dx_tbl_treated), file = "table_dx_num.tex", include.rownames = FALSE, include.colnames = TRUE,
      floating = TRUE, type = "latex", table.placement = "ht",
      floating.environment = "table",
      caption.placement = "top",
      caption = "Descriptive Information about Participants",
      label = "tab:dx")



################################################
######### Figure 6: Point Range Plots ##########
################################################

## Isolate treated subjects
post <- subset(df, Treated==1)

## Calculate lm() and arrange
pointranges <- post %>% group_by(T_Newtech,T_Specific,T_Political) %>% 
  do(data.frame(tidy(lm(Y_Likert ~ T_Specific*T_Newtech + T_Political, data = .),conf.int=T, conf.level = 0.95 ))) # 0.95 in paper which corresponds to p=0.001; should be 0.83 for p=0.05
pointranges <- subset(pointranges, !is.na(estimate))

## Factorize treatment conditions
pointranges$treatment <- ifelse(pointranges$T_Newtech==1, "new", "old")
pointranges$treatment <- ifelse(pointranges$T_Specific==1, paste0("spec", sep="-", pointranges$treatment), paste0("gen", sep="-", pointranges$treatment))
pointranges$treatment <- factor(pointranges$treatment, levels = c("spec-old", "gen-new",  "spec-new", "gen-old"), labels = c("spec-old", "gen-old",  "spec-new", "gen-new"))
pointranges$T_Political <- factor(pointranges$T_Political, labels = c("Neutral","Political"))

## Adding facet labels for plot
facet_labels <- c("spec-old" = "Mundane\nSpecific",
                  "gen-old" = "Mundane\nGeneral",
                  "spec-new" = "Novel\nSpecific",
                  "gen-new" = "Novel\nGeneral")

## Generate synthetic comparisons next to actual estimated data 
## **** (DO NOT CONFUSE SYNTHETIC DATA FOR ACTUAL ESTIMATES) ****
summarize <- dplyr::summarise 

pointranges_aug <- pointranges %>%
  group_by(treatment, T_Political) %>%
  summarize(
    estimate = mean(estimate), 
    conf.low = mean(conf.low), 
    conf.high = mean(conf.high), 
    .groups = 'drop'
  ) %>%
  ungroup() %>%
  # Duplicate the data for two types of fake points and create synthetic points
  bind_rows(
    pointranges %>%
      group_by(treatment, T_Political) %>%
      summarize(estimate = mean(estimate), conf.low = mean(conf.low), conf.high = mean(conf.high), .groups = 'drop') %>%
      mutate(T_Political = case_when(
        T_Political == "Political" ~ "Political (Priors)",
        T_Political == "Neutral" ~ "Neutral (Priors)",
        TRUE ~ as.character(T_Political)
      )),
    pointranges %>%
      group_by(treatment, T_Political) %>%
      summarize(estimate = mean(estimate) * 0.9, conf.low = mean(conf.low) * 0.9, conf.high = mean(conf.high) * 0.9, .groups = 'drop') %>%
      mutate(T_Political = case_when(
        T_Political == "Political" ~ "Political (Theory)",
        T_Political == "Neutral" ~ "Neutral (Theory)",
        TRUE ~ as.character(T_Political)
      ))
  )

## Add null and hypo point estimates and confidence intervals on synthetic data
pointranges_aug$estimate <- ifelse(grepl("spec-old", pointranges_aug$treatment) & pointranges_aug$T_Political=="Political (Priors)", 2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("spec-old", pointranges_aug$treatment) & pointranges_aug$T_Political=="Neutral (Priors)", 2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("gen-old", pointranges_aug$treatment) & pointranges_aug$T_Political=="Political (Priors)", -2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("gen-old", pointranges_aug$treatment) & pointranges_aug$T_Political=="Neutral (Priors)", 0, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("spec-new", pointranges_aug$treatment) & pointranges_aug$T_Political=="Political (Priors)", 2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("spec-new", pointranges_aug$treatment) & pointranges_aug$T_Political=="Neutral (Priors)", 2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("gen-new", pointranges_aug$treatment) & pointranges_aug$T_Political=="Political (Priors)", -2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("gen-new", pointranges_aug$treatment) & pointranges_aug$T_Political=="Neutral (Priors)", -2.1, pointranges_aug$estimate)

pointranges_aug$estimate <- ifelse(grepl("spec-old", pointranges_aug$treatment) & pointranges_aug$T_Political=="Political (Theory)", 2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("spec-old", pointranges_aug$treatment) & pointranges_aug$T_Political=="Neutral (Theory)", 2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("gen-old", pointranges_aug$treatment) & pointranges_aug$T_Political=="Political (Theory)", 0.5, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("gen-old", pointranges_aug$treatment) & pointranges_aug$T_Political=="Neutral (Theory)", 0.6, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("spec-new", pointranges_aug$treatment) & pointranges_aug$T_Political=="Political (Theory)", -2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("spec-new", pointranges_aug$treatment) & pointranges_aug$T_Political=="Neutral (Theory)", -2.1, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("gen-new", pointranges_aug$treatment) & pointranges_aug$T_Political=="Political (Theory)", 0.5, pointranges_aug$estimate)
pointranges_aug$estimate <- ifelse(grepl("gen-new", pointranges_aug$treatment) & pointranges_aug$T_Political=="Neutral (Theory)", 0.6, pointranges_aug$estimate)

pointranges_aug$conf.low <- ifelse(grepl("Priors|Theory", pointranges_aug$T_Political), pointranges_aug$estimate+.01, pointranges_aug$conf.low)
pointranges_aug$conf.high <- ifelse(grepl("Priors|Theory", pointranges_aug$T_Political), pointranges_aug$estimate-.01, pointranges_aug$conf.high)


## Update factor levels to include new synthetic points
pointranges_aug$T_Political <- factor(pointranges_aug$T_Political, levels = c(
  "Neutral (Priors)", "Political (Priors)", "Neutral (Theory)", "Political (Theory)", "Neutral", "Political"
))

## T3 on the legend
pointranges_aug <- pointranges_aug %>%
  mutate(
    T_Political = forcats::fct_recode(T_Political,
                                      "Political (ATE)" = "Political",
                                      "Neutral (ATE)" = "Neutral",
                                      "Political (Priors)" = "Political (Priors)",
                                      "Neutral (Priors)" = "Neutral (Priors)",
                                      "Political (Theory)" = "Political (Theory)",
                                      "Neutral (Theory)" = "Neutral (Theory)")
  )

## Custom shapes/colors for distinguishing *actual* from synthetic results
shape_settings <- c(
  "Political (Priors)" = 25, # Solid triangle for conventional wisdom
  "Neutral (Priors)" = 25,   # Solid triangle for conventional wisdom
  "Political (Theory)" = 24, # Solid square for theory
  "Neutral (Theory)" = 24,    # Solid square for theory
  "Political (ATE)" = 16,       # Solid circle for actual estimates
  "Neutral (ATE)" = 16         # Solid circle for actual estimates
)

color_settings <- c(
  "Political (ATE)" = "black",   
  "Neutral (ATE)" = "gray",        
  "Political (Priors)" = "black", 
  "Neutral (Priors)" = "gray",  
  "Political (Theory)" = "black", 
  "Neutral (Theory)" = "gray"    
)

pointranges_aug <- pointranges_aug %>%
  mutate(group = case_when(
    grepl("(Priors)", T_Political) ~ "Priors",
    grepl("(Theory)", T_Political) ~ "Theory",
    TRUE ~ "ATE" 
  ))

pointranges_aug <- pointranges_aug %>%
  mutate(group = factor(group, levels = c("Priors", "Theory", "ATE")))

facet_labels <- c("spec-old" = "Mundane Technology\n× Specific Law",
                  "gen-old" = "Mundane Technology\n× General Law",
                  "spec-new" = " Novel Technology\n× Specific Law",
                  "gen-new" = "Novel Technology\n× General Law")


## Plot synthetic data to the left of actual estimates
png(filename="ATE_pointrange.png", width=1200, height=300, res=160)

ggplot(pointranges_aug, aes(x = group, y = estimate, ymin = conf.low, ymax = conf.high,
                            shape = T_Political, color = T_Political)) +
  geom_pointrange(position = position_dodge(width = 0.5), size = 0.75) +
  facet_wrap(~treatment, nrow = 1, scales = "free_x", labeller = as_labeller(facet_labels)) +
  scale_shape_manual(values = shape_settings) +
  scale_color_manual(values = color_settings, breaks=c("Political (ATE)", "Neutral (ATE)")) +
  scale_x_discrete(labels = c("Priors", "Theory", "ATE")) + 
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(colour = "black"),
    panel.border = element_rect(colour = "black", fill = NA, size = 0.8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  ) +
  guides(
    color = guide_legend(title = "T3: Client Type"),
    shape = "none" 
  )

dev.off() 
graphics.off() # reset margins




################################################
##### Figure A.2: Randomization Inference ######
################################################

## Create an interaction variable
ridf <- subset(df, Treated==1) %>%
  mutate(Interaction = interaction(T_Newtech, T_Specific))

## Compute observed statistics for main effects
observed_model <- lm(Y_Likert ~ T_Newtech * T_Specific + T_Political, data = ridf)

## Extract coefficients for the main effects and interaction
observed_stat_newtech <- coef(observed_model)["T_Newtech"]
observed_stat_specific <- coef(observed_model)["T_Specific"]
observed_stat_political <- coef(observed_model)["T_Political"]
observed_stat_interaction <- coef(observed_model)["T_Newtech:T_Specific"]

## Function to permute treatments and compute test statistics
randomize <- function(data) {
  # Permute the treatment assignments
  permuted_data <- data %>%
    mutate(
      T_Newtech = sample(T_Newtech),
      T_Specific = sample(T_Specific),
      T_Political = sample(T_Political)  # Permute without interactions
    )
  
  # Fit the model with permuted data
  permuted_model <- lm(Y_Likert ~ T_Newtech * T_Specific + T_Political, data = permuted_data)
  
  # Extract coefficients for main effects and interaction
  permuted_stat_newtech <- coef(permuted_model)["T_Newtech"]
  permuted_stat_specific <- coef(permuted_model)["T_Specific"]
  permuted_stat_political <- coef(permuted_model)["T_Political"]
  permuted_stat_interaction <- coef(permuted_model)["T_Newtech:T_Specific"]
  
  return(c(permuted_stat_newtech, permuted_stat_specific, permuted_stat_political, permuted_stat_interaction))
}

## Generate null distributions for the test statistics
set.seed(888)
n_permutations <- 50000
null_distribution <- replicate(n_permutations, randomize(ridf))

## Extract null distributions for each effect
null_newtech <- null_distribution[1, ]
null_specific <- null_distribution[2, ]
null_political <- null_distribution[3, ]
null_interaction <- null_distribution[4, ]

## Calculate p-values for each effect
p_value_newtech <- mean(abs(null_newtech) >= abs(observed_stat_newtech))
p_value_specific <- mean(abs(null_specific) >= abs(observed_stat_specific))
p_value_political <- mean(abs(null_political) >= abs(observed_stat_political))
p_value_interaction <- mean(abs(null_interaction) >= abs(observed_stat_interaction))

## Print results
cat("Observed Statistic (T_newtech):", observed_stat_newtech, "\n")
cat("P-value (T_newtech):", p_value_newtech, "\n")

cat("Observed Statistic (T_specific):", observed_stat_specific, "\n")
cat("P-value (T_specific):", p_value_specific, "\n")

cat("Observed Statistic (T_political):", observed_stat_political, "\n")
cat("P-value (T_political):", p_value_political, "\n")

cat("Observed Statistic (Interaction):", observed_stat_interaction, "\n")
cat("P-value (Interaction):", p_value_interaction, "\n")


## Visualization function
plot_null <- function(null_distribution, observed_stat, title, xlab) {
  ggplot(data = data.frame(null_stat = null_distribution), aes(x = null_stat)) +
    geom_histogram(color = "black", fill = "lightblue", bins = 30) +
    geom_vline(xintercept = observed_stat, color = "red", linetype = "dashed", size = 1) +
    labs(title = title, x = xlab, y = "Frequency") +
    theme_minimal()
}

## Create plots for each null distribution
plot_newtech <- plot_null(null_newtech, observed_stat_newtech, "Null Distribution for T1: Novel", "Coefficient")
plot_specific <- plot_null(null_specific, observed_stat_specific, "Null Distribution for T2: Specific", "Coefficient")
plot_political <- plot_null(null_political, observed_stat_political, "Null Distribution for T3: Political", "Coefficient")
plot_interaction <- plot_null(null_interaction, observed_stat_interaction, "Null Distribution for T1xT2 Interaction", "Coefficient")


## Arrange plots in a grid
library(cowplot)

png(filename="ATE_RandomizationInference.png", width=1200, height=600, res=120)
plot_grid(
  plot_newtech, plot_specific, plot_political, plot_interaction,
  ncol = 2,  # Number of columns
  labels = c("A", "B", "C", "D"),  # Optional: Label each plot
  label_size = 12  # Optional: Set label size
)
dev.off()




################################################
########## Figure A.3: Balance Tests ########### 
################################################

library(cobalt)
library(WeightIt)

class(df)
balancedf <- data.frame(sapply(df, function(x) as.numeric(as.character(x))))
sapply(balancedf, class)

## Exclude NAs
balancedf <- subset(balancedf, !is.na(Y_Likert)) # only if treated - should already be selected on this
balancedf <- balancedf %>% dplyr::select(where(~ !all(is.na(.)))) # delete cols that only contain NAs (numeric conversion)

## Exclude aversion variables (unavailable for most of sample)
balancedf <- balancedf %>% dplyr::select(-matches("Aversion"))
x_columns <- grep("^(T|X|Y)_|Wave", names(balancedf), value = TRUE)
balancedf <- balancedf[, x_columns]

## Confirm no NAs
colnames(balancedf)[colSums(is.na(balancedf)) > 0]

## Check again for missingness
balancedf %>%
  dplyr::summarise(across(everything(), ~ any(is.na(.)))) %>%
  pivot_longer(everything(), names_to = "column", values_to = "has_missing") %>%
  filter(has_missing) %>%
  pull(column)

## Remove columns that contain any NA values
balancedf <- balancedf[, colSums(is.na(balancedf)) == 0]

## Create main interaction term
names(balancedf)
balancedf$T1xT2 <- balancedf$T_Newtech + balancedf$T_Specific + (balancedf$T_Newtech*balancedf$T_Specific)


names(balancedf)
treats <- balancedf %>% dplyr::select(matches("T_|T1xT2"))
outs <- balancedf %>% dplyr::select(matches("Y_Likert"))
covs <- balancedf %>% dplyr::select(matches("X_"))

## Print results
b1 <- bal.tab(x = covs, treat = treats$T_Newtech, int = TRUE)
b2 <- bal.tab(x = covs, treat = treats$T_Specific, int = TRUE)
b3 <- bal.tab(x = covs, treat = treats$T_Political, int = TRUE)

## Create a factor variable from combinations of col1, col2, col3
names(covs)
treats$T_T1xT2 <- interaction(treats$T_Newtech, treats$T_Specific, sep = "-")
b4  <- bal.tab(x = covs, treat = treats$T_T1xT2, int = TRUE)

## Assign weights
W.new <- weightit(T_Newtech ~ covs, data = treats, method = "ps", estimand = "ATE", maxit=100)
W.spec <- weightit(T_Specific ~ covs, data = treats, method = "ps", estimand = "ATE", maxit=100) 
W.pol <- weightit(T_Political ~ covs, data = treats, method = "ps", estimand = "ATE", maxit=100)
W.t1t2 <- weightit(T_T1xT2 ~ covs, data = treats, method = "ps", estimand = "ATE", maxit=100, use.mlogit = F)

v1 <- var.names(b1, type = "vec", minimal = TRUE)
v2 <- var.names(b2, type = "vec", minimal = TRUE)
v3 <- var.names(b3, type = "vec", minimal = TRUE)
v4 <- var.names(b4, type = "vec", minimal = TRUE)

bal.tab(W.spec)
bal.tab(W.new)
bal.tab(W.pol)
bal.tab(W.t1t2)


## T3
l1 <- love.plot(W.new, stars="raw", colors=c("black","gray50"), 
                var.names = v1, threshold = c(m = .25, ks = .05),  # m CI default = .1
                limits = c(-.4, 1), grid = FALSE, size=4,
                shapes = c("triangle filled", "circle filled"),
                line = TRUE, drop.distance = T, alpha=.8,
                binary = "std", title = "T1: Novel",
                drop.missing = T,
                position = c(.75, .25)) +
  theme(legend.box.background = element_rect())

## T2
l2 <- love.plot(W.spec, stars="raw", colors=c("black","gray50"), 
                var.names = v2,  threshold = c(m = .25, ks = .05),  # m CI default = .1
                limits = c(-.4, 1), grid = FALSE, size=4,
                shapes = c("triangle filled", "circle filled"),
                line = TRUE, drop.distance = T, alpha=.8,
                binary = "std", title = "T2: Specific",
                drop.missing = T,
                position = c(.75, .25)) +
  theme(legend.box.background = element_rect())


## T3
l3 <- love.plot(W.pol, stars="raw", colors=c("black","gray50"), 
                var.names = v3,  threshold = c(m = .25, ks = .05),  # m CI default = .1
                limits = c(-.4, 1), grid = FALSE, size=4,
                shapes = c("triangle filled", "circle filled"),
                line = TRUE, drop.distance = T, alpha=.8,
                binary = "std", title = "T3: Political Client",
                position = c(.75, .25)) +
  theme(legend.box.background = element_rect())

## T1 x T2
l4 <- love.plot(W.t1t2, stars="raw", colors=c("black","gray50"), 
                var.names = v4,  threshold = c(m = .25, ks = .05),  # m CI default = .1
                limits = c(-.4, 1), grid = FALSE, size=4,
                shapes = c("triangle filled", "circle filled"),
                line = TRUE, drop.distance = T, alpha=.8,
                binary = "std", title = "T1 x T2",
                position = c(.75, .25)) +
  theme(legend.box.background = element_rect())


## Plot balance on observables
png(filename="ATE_BalancePlots.png", width=1100, height=1500, res=150)
plot_grid(l1, l2, l3, l4,
          nrow = 2,
          labels = "AUTO",
          label_size = 12,
          align = "v"
)

dev.off() 
graphics.off() # reset margins




################################################
######### Figure A.4: LASSO CV/Boostrap ######## 
################################################

library(glmnet)
library(caret)
library(dplyr)

set.seed(888)

## Train/test split (90/10)
split <- createDataPartition(y = balancedf$Y_Likert, p = 0.9, list = FALSE)

X_train <- balancedf[split, ] %>%
  dplyr::select(-matches("Y_|X_LSAT_Score|Wave")) %>%
  as.matrix()

Y_train <- as.numeric(as.character(balancedf$Y_Likert[split]))

X_test <- balancedf[-split, ] %>%
  dplyr::select(-matches("Y_|X_LSAT_Score|Wave")) %>%
  as.matrix()

Y_test <- as.numeric(as.character(balancedf$Y_Likert[-split]))

cat("Length of Y_test:", length(Y_test), "\n")
cat("Length of X_test:", nrow(X_test), "\n")

## Elastic-net grid (alpha from 0 to 1)
alpha_values <- seq(0, 1, length.out = 8)  # 0=ridge, 1=lasso

## Fixed CV folds (deterministic)
set.seed(888)
foldid <- sample(rep(1:10, length.out = length(Y_train)))

## Fixed lambda path (deterministic lambda)
lambda_grid <- exp(seq(log(1e-4), log(100), length.out = 200)) # adjust range if desired

## Fit cv.glmnet for each alpha
cv_results <- lapply(alpha_values, function(a) {
  cv.glmnet(
    x = X_train,
    y = Y_train,
    alpha = a,
    foldid = foldid,
    lambda = lambda_grid,
    standardize = TRUE
  )
})
names(cv_results) <- paste0("cv_", alpha_values)

## Pick best alpha by min CV MSE
min_errors <- sapply(cv_results, function(m) min(m$cvm))
best_alpha <- alpha_values[which.min(min_errors)]
best_model <- cv_results[[paste0("cv_", best_alpha)]]

best_lambda <- best_model$lambda.min  # or best_model$lambda.1se
cat("Best alpha:", best_alpha, "\n")
cat("Minimum CV MSE:", min(min_errors), "\n")
cat("Best lambda:", best_lambda, "\n")


## Extract coefficients at best lambda
best_coefs <- coef(best_model, s = best_lambda)
print(best_coefs)


## Generate OOS predictions using the best model
predictions <- predict(best_model, s = best_lambda, newx = X_test)

## Calculate Mean Squared Error (MSE) on the test set
mse <- mean((Y_test - predictions)^2)

## Larger test MSE implies overfitting, smaller implies underfitting.
cat("MSE on training data (min CV MSE):", min(min_errors), "\n")  # LOW => good in-sample predictive performance
cat("MSE on test data:", mse, "\n")                               # CLOSE => no overfitting


## Bootstrapped CIs
library(boot)

n_boot <- 50000

dummy_data <- seq_len(nrow(X_train))

boot_func <- function(data, indices) {
  
  ## Lock down RNG per replicate
  set.seed(888 + indices[1])
  
  n <- nrow(X_train)
  
  ## Bootstrap resample for training
  train_indices <- sample.int(n, size = round(0.7 * n), replace = TRUE)
  if (length(train_indices) < 10) return(NA_real_)
  
  test_indices <- setdiff(seq_len(n), train_indices)
  if (length(test_indices) < 2) return(NA_real_)
  
  X_train_boot <- X_train[train_indices, , drop = FALSE]
  Y_train_boot <- Y_train[train_indices]
  X_test_boot  <- X_train[test_indices, , drop = FALSE]
  Y_test_boot  <- Y_train[test_indices]
  
  if (var(Y_train_boot, na.rm = TRUE) == 0) return(NA_real_)
  
  ## Fit model with fixed alpha & lambda
  model_fold <- glmnet(
    x = X_train_boot,
    y = Y_train_boot,
    alpha = best_alpha,
    lambda = best_lambda,
    standardize = TRUE
  )
  
  ## Test MSE for this bootstrap sample
  test_pred <- predict(model_fold, newx = X_test_boot, s = best_lambda)
  boot_test_mse <- mean((Y_test_boot - test_pred)^2)
  
  ## Return Test MSE − CV Error
  boot_test_mse - min(min_errors)
}

set.seed(888)
boot_results <- boot(
  data = dummy_data,
  statistic = boot_func,
  R = n_boot
)

## Remove NA reps
boot_vals <- boot_results$t
boot_vals <- boot_vals[!is.na(boot_vals)]

## Percentile CI
boot_ci <- quantile(boot_vals, probs = c(0.025, 0.975), names = FALSE)

cat(
  "Bootstrap 95% CI for Test MSE − CV Error:",
  boot_ci, "\n"
)

## ============================================================
## LASSO / Elastic-Net Cross-Validation Plots
## ============================================================

png(filename = "Lasso_MSE.png", width = 1000, height = 800, res = 120)

par(mfrow = c(3, 3), mar = c(1, 3, 7, 2))

## Individual CV curves (one per alpha)
for (i in seq_along(alpha_values)) {
  alpha <- alpha_values[i]
  model <- cv_results[[paste0("cv_", alpha)]]
  plot(model, main = paste("alpha =", round(alpha, 3)))
}

dev.off()


png(filename = "Lasso_MSE.png", width = 1000, height = 800, res = 120)
par(mfrow = c(3, 3), mar = c(1, 3, 7, 2))

## 8 individual CV plots
for (i in seq_along(alpha_values)) {
  alpha <- alpha_values[i]
  model <- cv_results[[paste0("cv_", alpha)]]
  plot(model, main = paste("alpha =", round(alpha, 3)))
}

## Overlay in the 9th panel
line_types  <- c("dotted", "dotdash", "dashed", "twodash",
                 "longdash", "solid", "solid", "solid")
line_widths <- c(1, 1, 1, 1, 1, 1, 1, 2)

plot(
  log(cv_results[[1]]$lambda),
  cv_results[[1]]$cvm,
  type = "l",
  lty = line_types[1],
  lwd = line_widths[1],
  xlim = c(-3, 7),
  ylim = range(sapply(cv_results, function(m) m$cvm)),
  xlab = "log(Lambda)",
  ylab = cv_results[[1]]$name,
  main = "CV MSE Overlay Across α"
)

for (i in 2:length(alpha_values)) {
  model <- cv_results[[paste0("cv_", alpha_values[i])]]
  lines(
    log(model$lambda),
    model$cvm,
    lty = line_types[i],
    lwd = line_widths[i]
  )
}

legend(
  "topright",
  legend = paste("alpha =", round(alpha_values, 3)),
  cex = 0.85,
  bty = "n",
  lty = line_types,
  lwd = line_widths
)
dev.off()





################################################
######### Table A.1: Numerical Outcomes ######## 
################################################

## 1. Fit base model
summary(fit.base <- lm(Y_Likert~T_Specific*T_Newtech+T_Political, data=balancedf)) 

## 2. Fit 3W interaction model
summary(fit.base.int <- lm(Y_Likert~T_Specific*T_Newtech*T_Political, data=balancedf)) 

## 3. Extract Lasso results 
exclude_vars <- c("(Intercept)", "T_Specific", "T_Newtech", "T_Political", "T1xT2")
tol <- 1e-6  # helps avoid tiny numerical jitter changing selection

significant_vars <- rownames(best_coefs)[
  abs(best_coefs[, 1]) > tol &
    !is.na(best_coefs[, 1]) &
    !(rownames(best_coefs) %in% exclude_vars)
]

cat("Selected vars (post-LASSO):\n")
print(significant_vars)

# Handle case where LASSO selects additional vars
rhs_extra <- if (length(significant_vars) > 0) {
  paste(significant_vars, collapse = " + ")
} else {
  "1"  # just keep intercept; base terms added below 
}

formula <- as.formula(
  paste("Y_Likert ~ T_Specific * T_Newtech + T_Political +", rhs_extra)
)

## 3. Fit the LASSO model
summary(fit.lasso <- lm(formula, data = balancedf))

## 4. Fit a subgroup model on covariates of special theoretical interest, searching for politicization in the "ambiguous" condition
summary(fit.evade <- lm(Y_Likert ~ T_Specific * T_Newtech + T_Political*(X_CNPractice + X_Polity + X_Conservatism_N + X_Isolationism_N + X_Practitioner), data = balancedf))
summary(fit.evade <- lm(Y_Likert ~ T_Newtech + T_Political*(X_From_Autocracy + X_Conservatism_N + X_Isolationism_N + X_Practitioner), data = subset(balancedf, T_Specific==0)))


## Plot numerical results
library(stargazer)
stargazer(fit.base, fit.base.int, fit.lasso, fit.evade,
          type = "latex",
          font.size = "small",
          column.sep.width = "-5pt",
          align = TRUE,
          column.labels = c("OLS-Unadjusted", "OLS 3W Interaction", "LASSO-Adjusted", "General Only"),
          single.row = TRUE,
          no.space = TRUE,
          object.names = FALSE,
          model.numbers = TRUE,
          dep.var.labels = c("Perceived Restrictiveness of Article II"),
          intercept.top = TRUE,
          intercept.bottom = FALSE,
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes.append = FALSE,
          notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"))




################################################
######## Figure A.5: Secondary Outcomes ######## 
################################################

library(ggpattern)

## Confidence ~ Specificity
p1 <- ggplot(post, aes(x = Y_Confidence, fill = T_Specific)) +
  # Hatching for level "1"
  geom_density_pattern(
    data = subset(post, T_Specific == "1"),
    aes(pattern = T_Specific),
    pattern = "stripe", # stripe, ellipse, crosshatch
    fill = "black", # Fill color for stripes
    colour = "black", # Border color for stripes
    pattern_density = 0.01, # Density of stripes
    pattern_spacing = 0.05, # Spacing between stripes
    pattern_angle = 135, # Angle of stripes
    alpha = 0.5 # Transparency of the pattern fill
  ) +
  # Solid fill for level "0"
  geom_density(
    data = subset(post, T_Specific == "0"),
    alpha = 0.5,
    fill = "gray"
  ) +
  scale_fill_manual(values = c("0" = "gray", "1" = "black")) +
  theme_minimal() +
  labs(
    fill = "Specificity",
    x = "Confidence",
    y = "Density",
    title = "Confidence Density",
    subtitle = ""
  )

## Confidence ~ Specificity
p2 <- ggplot(post, aes(x = Y_CitePreamble, fill = T_Specific)) +
  # Hatching for level "1"
  geom_density_pattern(
    data = subset(df, Treated==1 & T_Specific == "1"),
    aes(pattern = T_Specific),
    pattern = "stripe", # stripe, ellipse, crosshatch
    fill = "black", # Fill color for stripes
    colour = "black", # Border color for stripes
    pattern_density = 0.01, # Density of stripes
    pattern_spacing = 0.05, # Spacing between stripes
    pattern_angle = 135, # Angle of stripes
    alpha = 0.5 # Transparency of the pattern fill
  ) +
  # Solid fill for level "0"
  geom_density(
    data = subset(df, Treated==1 & T_Specific == "0"),
    alpha = 0.5,
    fill = "gray"
  ) +
  scale_fill_manual(values = c("0" = "gray", "1" = "black")) +
  theme_minimal() +
  labs(
    fill = "Specificity",
    x = "Attention to Spirit",
    y = "Density",
    title = "Attention to Spirit",
    subtitle = ""
  )

## Confidence ~ Specificity
p3 <- ggplot(subset(df, Treated==1), aes(x = Y_CiteDefinitions, fill = T_Specific)) +
  # Hatching for level "1"
  geom_density_pattern(
    data = subset(df, Treated==1 & T_Specific == "1"),
    aes(pattern = T_Specific),
    pattern = "stripe", # stripe, ellipse, crosshatch
    fill = "black", # Fill color for stripes
    colour = "black", # Border color for stripes
    pattern_density = 0.01, # Density of stripes
    pattern_spacing = 0.05, # Spacing between stripes
    pattern_angle = 135, # Angle of stripes
    alpha = 0.5 # Transparency of the pattern fill
  ) +
  # Solid fill for level "0"
  geom_density(
    data = subset(df, Treated==1 & T_Specific == "0"),
    alpha = 0.5,
    fill = "gray"
  ) +
  scale_fill_manual(values = c("0" = "gray", "1" = "black")) +
  theme_minimal() +
  labs(
    fill = "Specificity",
    x = "Attention to Definitions",
    y = "Density",
    title = "Attention to Definitions",
    subtitle = ""
  )

## Wordcount ~ Specificity
p4 <- ggplot(post, aes(x = log(Y_Nchar), fill = T_Specific)) +
  # Hatching for level "1"
  geom_density_pattern(
    data = subset(post, T_Specific == "1"),
    aes(pattern = T_Specific),
    pattern = "stripe", # stripe, ellipse, crosshatch
    fill = "black", # Fill color for stripes
    colour = "black", # Border color for stripes
    pattern_density = 0.01, # Density of stripes
    pattern_spacing = 0.05, # Spacing between stripes
    pattern_angle = 135, # Angle of stripes
    alpha = 0.5 # Transparency of the pattern fill
  ) +
  # Solid fill for level "0"
  geom_density(
    data = subset(post, T_Specific == "0"),
    alpha = 0.5,
    fill = "gray"
  ) +
  scale_fill_manual(values = c("0" = "gray", "1" = "black")) +
  theme_minimal() +
  labs(
    fill = "Specificity",
    x = "Wordcount",
    y = "Density",
    title = "Log Wordcount Density",
    subtitle = ""
  )

## Edits ~ Specificity
p5 <- ggplot(post, aes(x = log(Timer_Clicks), fill = T_Specific)) +
  # Hatching for level "1"
  geom_density_pattern(
    data = subset(post, T_Specific == "1"),
    aes(pattern = T_Specific),
    pattern = "stripe", # stripe, ellipse, crosshatch
    fill = "black", # Fill color for stripes
    colour = "black", # Border color for stripes
    pattern_density = 0.01, # Density of stripes
    pattern_spacing = 0.05, # Spacing between stripes
    pattern_angle = 135, # Angle of stripes
    alpha = 0.5 # Transparency of the pattern fill
  ) +
  # Solid fill for level "0"
  geom_density(
    data = subset(post, T_Specific == "0"),
    alpha = 0.5,
    fill = "gray"
  ) +
  scale_fill_manual(values = c("0" = "gray", "1" = "black")) +
  theme_minimal() +
  labs(
    fill = "Specificity",
    x = "Log # Edits",
    y = "Density",
    title = "Log Number of Edits",
    subtitle = "'"
  )

## Time Spent ~ Specificity (controlling for typing speed)
p6 <- ggplot(post, aes(x = log(Timer_Write), fill = T_Specific)) +
  # Hatching for level "1"
  geom_density_pattern(
    data = subset(post, T_Specific == "1"),
    aes(pattern = T_Specific),
    pattern = "stripe", # stripe, ellipse, crosshatch
    fill = "black", # Fill color for stripes
    colour = "black", # Border color for stripes
    pattern_density = 0.01, # Density of stripes
    pattern_spacing = 0.05, # Spacing between stripes
    pattern_angle = 135, # Angle of stripes
    alpha = 0.5 # Transparency of the pattern fill
  ) +
  # Solid fill for level "0"
  geom_density(
    data = subset(post, T_Specific == "0"),
    alpha = 0.5,
    fill = "gray"
  ) +
  scale_fill_manual(values = c("0" = "gray", "1" = "black")) +
  theme_minimal() +
  labs(
    fill = "Specificity",
    x = "Log Time Spent",
    y = "Density",
    title = "Log Time Spent",
    subtitle = "'"
  )



## Plot secondary outcomes
library(cowplot)

## X ~ as Rule Specificity Varies
## Visualized using hatching for '1' and solid fill for '0
png(filename="ATE_OtherOutcomes.png", width=1200, height=600, res=120)
plot_grid(p1, p2, p3, p4, p5, p6,
          nrow = 3,
          labels = "AUTO",
          label_size = 12,
          align = "v"
)
dev.off()

## Relationship with T2_Specific, accounting for LSAT Score/Cumulative Years of Training
## Note: X_Edu_Completed_N is capped at 1 year postgrad
summary(fit <- lm(Y_Nchar ~ T_Specific * (X_LSAT_Score + X_Edu_Completed_N), data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))
summary(fit <- lm(Y_Confidence ~ T_Specific * (X_LSAT_Score + X_Edu_Completed_N), data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))




################################################
######### Figure A.6: Sample Geography ######### 
################################################

library(viridis)
library(countrycode)
library(rnaturalearth)

# Load world map, exclude Antarctica
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>%
  filter(name != "Antarctica")

df$X_LocationISO <- countrycode(df$X_Location_From, origin = "country.name", destination = "iso3c")

# Countries represented in dataset
country_counts <- df %>%
  dplyr::group_by(X_LocationISO) %>%
  dplyr::summarise(N_participants = n())

map_data <- world %>%
  left_join(country_counts, by = c("iso_a3" = "X_LocationISO"))

## Plot outbound geographical composition
png(filename="samplemap.png", width=2400, height=1200, res=300)
ggplot(map_data) +
  geom_sf(aes(fill = log(N_participants))) +
  scale_fill_gradient(low = "lightblue", high = "navyblue",
                      name = "Log #", na.value = "grey90") +
  theme_void() +
  labs(title = "") +
  theme(
    legend.position = "right",
    legend.position.inside = c(0.85, 0.5), 
    plot.title = element_text(size = 20),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 16)
  )
dev.off()

# Remove helper variable
df$X_LocationISO <- NULL

# 64 countries represented
length(unique(df$X_Location_From))

# Range of “lawyerliness” of sampled countries 
range(na.omit(df$X_LawPtP))





################################################
######### Figure A.7: Wave Stability  ##########
################################################

library(dplyr)
library(broom)
library(ggplot2)
library(stringr)
library(car)

# Sort schools alphabetically or numerically and assign an order
df <- df[order(df$EndDate), ]
df$WaveOrder <- seq_len(nrow(df))

## Define baseline model function
run_model <- function(data) {
  lm(Y_Likert ~ T_Newtech*T_Specific + T_Political + WaveOrder + T_Newtech*T_Specific,
     data = data)
}

## Get list of waves
waves <- unique(df$Wave)

## Run leave-one-wave-out models
jackknife_results <- lapply(waves, function(s) {
  model <- run_model(subset(df, Wave != s))
  coefs <- broom::tidy(model) %>%
    mutate(left_out = s)
  return(coefs)
}) %>%
  bind_rows()


## Filter for the key coefficient
plot_data <- jackknife_results %>%
  filter(term == "T_Newtech:T_Specific") %>%
  mutate(left_out = factor(left_out, levels = unique(left_out)))

## Estimate full-sample coefficient for reference line
full_model <- run_model(df)
full_coef <- coef(full_model)["T_Newtech:T_Specific"]


## Plot
png(filename="ATE_jackknife.png", width=1200, height=400, res=140)
ggplot(plot_data, aes(x = left_out, y = estimate)) +
  geom_point(size = 3, color = "darkblue") +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
                width = 0.15, color = "steelblue") +
  geom_hline(yintercept = full_coef, linetype = "dashed", color = "red", linewidth = 0.8) +
  labs(
    title = "",
    subtitle = "Coefficient for T_Newtech × T_Specific",
    x = "Omitted Wave",
    y = "Estimate (±95% CI)"
  ) +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))
dev.off()

# Jackknife variance ratio
jack_var <- var(jackknife_results$estimate[jackknife_results$term == "T_Newtech:T_Specific"], na.rm = TRUE)
mean_se2 <- mean(jackknife_results$std.error[jackknife_results$term == "T_Newtech:T_Specific"]^2, na.rm = TRUE)
jack_ratio <- jack_var / mean_se2

# Print jackknife variance ratio
jack_ratio

## Delete helper column
df$WaveOrder <- NULL






################################################
###### Table A.3: Untreated Sample Stats #######
################################################

## Function to generate summaries for numeric and factor variables
dx_tbl_untreated <- summarize_data(subset(dx_tbl, Treated!=1 | is.na(Treated))[,1:51])

rownames(dx_tbl_untreated) <- NULL
xtable(dx_tbl_untreated) %>%
  print(include.rownames = FALSE)

## As percentages
round(prop.table(table(subset(dx_tbl_untreated, Treated==1)$LawSocID)) * 100, 2)

xtable(dx_tbl_untreated)

print(xtable(dx_tbl_untreated), file = "table_dx_untreated.tex", include.rownames = FALSE, include.colnames = TRUE,
      floating = TRUE, type = "latex", table.placement = "ht",
      floating.environment = "table",
      caption.placement = "top",
      caption = "Descriptive Information about Untreated Participants",
      label = "tab:dx")



################################################
#### Manipulation Checks for Treated Sample ####
################################################

df$Q_ManipTech_F
df$Q_ManipLaw_F
df$Q_ManipAware_F
df$Q_ManipInterp_F # not a great check - technically subjective

## Full sample
which(subset(df, Treated==1)$Q_ManipLaw_F==1 &
        subset(df, Treated==1)$Q_ManipTech_F==1 &
        subset(df, Treated==1)$Q_ManipAware_F==1)

## Without pilot
which(subset(df, Treated==1 & Wave!="PILOT 2020")$Q_ManipLaw_F==1 &
        subset(df, Treated==1 & Wave!="PILOT 2020")$Q_ManipTech_F==1 &
        subset(df, Treated==1 & Wave!="PILOT 2020")$Q_ManipAware_F==1)

## Each check individually
summary(subset(df, Treated==1)$Q_ManipLaw_F)
summary(subset(df, Treated==1)$Q_ManipTech_F)
summary(subset(df, Treated==1)$Q_ManipAware_F)
summary(subset(df, Treated==1)$Q_ManipInterp_F)

summary(subset(df, Treated==1 & Wave!="PILOT 2020")$Q_ManipLaw_F)
summary(subset(df, Treated==1 & Wave!="PILOT 2020")$Q_ManipTech_F)
summary(subset(df, Treated==1 & Wave!="PILOT 2020")$Q_ManipAware_F)
summary(subset(df, Treated==1 & Wave!="PILOT 2020")$Q_ManipInterp_F)

## As percentages
n <- nrow(subset(df, Treated == 1))

## Percent failing each check
p1 <- mean(subset(df, Treated == 1)$Q_ManipTech_F == 0, na.rm = TRUE) * 100
p2 <- mean(subset(df, Treated == 1)$Q_ManipLaw_F == 0, na.rm = TRUE) * 100
p3 <- mean(subset(df, Treated == 1)$Q_ManipAware_F == 0, na.rm = TRUE) * 100
p4 <- mean(subset(df, Treated == 1)$Q_ManipInterp_F == 0, na.rm = TRUE) * 100

cat(sprintf("Pass ManipCheck1: %.1f%%\n", p1))
cat(sprintf("Pass ManipCheck2: %.1f%%\n", p2))
cat(sprintf("Pass ManipCheck3: %.1f%%\n", p3))
cat(sprintf("Pass ManipCheck4: %.1f%%\n", p4))

## Joint failures
p_all <- mean(
  subset(df, Treated == 1)$Q_ManipTech_F == 1 &
    subset(df, Treated == 1)$Q_ManipLaw_F == 1 &
    subset(df, Treated == 1)$Q_ManipAware_F == 1,
  na.rm = TRUE
) * 100

cat(sprintf("Failed all objective checks: %.1f%%\n", p_all))

## Relationship between difficulty (proxied by aptitude) and pass rates
summary(fit <- lm(Q_ManipLaw_F ~ X_LSAT_Score * X_Edu_Completed_N, data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))
summary(fit <- lm(Q_ManipTech_F ~ X_LSAT_Score * X_Edu_Completed_N, data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))
summary(fit <- lm(Q_ManipInterp_F ~ X_LSAT_Score * X_Edu_Completed_N, data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))
summary(fit <- lm(Q_ManipAware_F ~ X_LSAT_Score * X_Edu_Completed_N, data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))

## Relationship between answer confidence and pass rates
summary(fit <- lm(Q_ManipLaw_F ~ Y_Confidence, data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))
summary(fit <- lm(Q_ManipTech_F ~ Y_Confidence, data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))
summary(fit <- lm(Q_ManipInterp_F ~ Y_Confidence, data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))
summary(fit <- lm(Q_ManipAware_F ~ Y_Confidence, data = subset(df, Treated==1 & !is.na(X_LSAT_Score))))




################################################
########## A.3: Heterogeneous Effects ########## 
################################################

rm(list=setdiff(ls(), "df"))

## Recommend reloading R 
loadedNamespaces()
sessionInfo()

## Java RAM settings https://www.oracle.com/java/technologies/javase/vmoptions-jsp.html
XLConnect::xlcFreeMemory() # frees java virtual memory
options(java.parameters = c("-XX:+UseConcMarkSweepGC", "-Xmx10g")) 
options(bartMachine_num_cores = parallel::detectCores() - 1) # Reduce to save memory
library(bartMachine)

## If doesn't take, launch from global options instead:
# 1. open terminal
# 2. type: nano ~/.zshrc
# 3. add line: export JAVA_TOOL_OPTIONS="-Xmx60g -XX:+UseG1GC"
# 4. add line open -a RStudio
# 5. click ctrl+o to save and ctrl+x to exit.
# 6. then, in terminal, run: source ~/.zshrc to launch rstudio w/ special options.

## Confirm working:
runtime_instance <- .jcall("java/lang/Runtime", "Ljava/lang/Runtime;", "getRuntime")
max_memory <- .jcall(runtime_instance, "J", "maxMemory") / 1024^3
print(paste("Max memory allocated to JVM:", max_memory, "GB"))
options("bartMachine_num_cores")

library(RcppParallel)
setThreadOptions(numThreads = parallel::detectCores() - 1)

library(dplyr)
library(tidyr)
library(stringr)
library(data.table)

## Clear memory
library(caret)
library(parallel)
set_bart_machine_num_cores(1) 
gc(verbose=TRUE, full=TRUE);gc(verbose=TRUE, full=TRUE);gc(verbose=TRUE, full=TRUE)
XLConnect::xlcFreeMemory()

df <- subset(df, !is.na(Y_Likert))

## Form BART df
X_full <- df %>% dplyr::select(matches("T1xT2|T_|X_|Y_Likert"))

## Remove factors (optional) 
unique(sapply(X_full, class))
names(X_full)[sapply(X_full, is.character)] 
names(X_full)[sapply(X_full, is.factor)]
X_full <- X_full[, sapply(X_full, class) != "factor"] # delete factors
X_full <- X_full[, sapply(X_full, class) != "character"] # delete factors

## Checking for NAs in the data before building the model
if(any(is.na(X_full)) || any(is.na(X_full))) {
  warning("NA values found in training data.")
} else {
  print("No NA values in training data. Proceeding with model build.")
}

names(X_full)[sapply(X_full, anyNA)]

## Get the names of all other columns excluding 'T_newtech' and 'T_specific'
treat_vars <- setdiff(names(X_full), c("T_Newtech", "T_Specific", "T_Political"))

## Limit to 2-way Treatment x Covariate interactions
## Loop through each of these other variables to create interaction terms
for(var in treat_vars) {
  # Create a new interaction column for each variable
  X_full[[paste("T1", var, sep = " x ")]] <- X_full$T_Newtech * X_full[[var]]
}

for(var in treat_vars) {
  # Create a new interaction column for each variable
  X_full[[paste("T2", var, sep = " x ")]] <- X_full$T_Specific * X_full[[var]]
}

for(var in treat_vars) {
  # Create a new interaction column for each variable
  X_full[[paste("T3", var, sep = " x ")]] <- X_full$T_Political * X_full[[var]]
}

X_full$`T1 x T2` <- X_full$T_Newtech * X_full$T_Specific

## Remove non-interactions
X_full_backup <- X_full
X_full <- X_full %>% dplyr::select(-matches("^(T|X)_"))

## Scale DV using scale() and adjust to range [-0.5, 0.5]
X_full$Y_Likert <- as.numeric(scale(X_full$Y_Likert, center = TRUE, scale = diff(c(-0.5, 0.5)) / 2)) # scale resp only

## Split into test/training (use 90/10 for better fit)
set.seed(888) 
df$randomizer <- sample.int(nrow(df), nrow(df), replace=F)

set.seed(888)
split <- createDataPartition(y = df$"randomizer", p=.9, list=FALSE) 

## Create response
X_train <- X_full[split,]; nrow(X_train)
Y_train <- as.numeric(as.character(X_train$Y_Likert))
X_train <- X_train %>% dplyr::select(matches("T_|X_|T1 x T2"))

X_test <- X_full[-split,]; nrow(X_test)
Y_test <- as.numeric(as.character(X_test$Y_Likert))
X_test <- X_test %>% dplyr::select(matches("T_|X_|T1 x T2"))


## Build a CV BART model
bart_machine_num_cores()
set_bart_machine_num_cores(1) 
gc(verbose=TRUE, full=TRUE);gc(verbose=TRUE, full=TRUE);gc(verbose=TRUE, full=TRUE)
XLConnect::xlcFreeMemory()

set_bart_machine_num_cores(11) 

treelist <- c(10,20,50,100,150,200) # or allow default

system.time({
  set.seed(888)
  bart_model_cv_man <- build_bart_machine_cv(as.data.frame(X_train), Y_train,
                                             use_missing_data = TRUE, 
                                             use_missing_data_dummies_as_covars = TRUE,
                                             mem_cache_for_speed = TRUE,
                                             flush_indices_to_save_RAM = TRUE,
                                             verbose = TRUE,
                                             # debug_log = FALSE,
                                             seed=888,
                                             serialize=TRUE,
                                             num_tree_cvs = treelist,
                                             k_cvs = c(2, 3, 5),
                                             k_folds = 20) # run m * K times
}); bart_model_cv_man # shapiro-wilk test for normality of the residuals

system.time({
  set.seed(888)
  var_importance_man <- investigate_var_importance(bart_model_cv_man, num_replicates_for_avg = 20)
}); var_importance_man


#saveRDS(bart_model_cv_man, "bart_model_cv_.rds")
#saveRDS(var_importance_man, file = "var_importance.rds")



######## DIAGNOSTICS ######## 

## Tree sensitivity 
set.seed(888)
tree_rmse <- rmse_by_num_trees(bart_model_cv_man, tree_list = treelist, plot = TRUE); tree_rmse # tree number

optimal_tree <- which.min(tree_rmse) # find lowest RMSE in output
optimal_tree <- treelist[optimal_tree] # use it to find optimal tree
# use this if want to rerun a non-cv model with optim_trees



################################################
##### Figure A.8: Interaction Importance #######
################################################

# Set the margin sizes (bottom, left, top, right)
graphics.off()
png(filename="BART_varimportance_manual.png", width=900, height=1200, res=90)
par(mar=c(4, 10, 4, 2) + 0.1)

# Order the Top 20 var_importance_props in descending order
var_importance_props <- var_importance_man$avg_var_props#[1:75]
se_var_importance <- var_importance_man$sd_var_props#[1:75]  # Replace with actual standard errors
ci_var_importance <- 1.96 * (se_var_importance / sqrt(nrow(X_train)))  # 95% confidence interval

# Set a scaling factor to avoid R treating very low SEs as zero-length
scaling_factor <- 1 

# Scale the confidence intervals
lower_bound <- var_importance_props - (ci_var_importance * scaling_factor)
upper_bound <- var_importance_props + (ci_var_importance * scaling_factor)

bp <- barplot(rev(var_importance_props),
              main = "Interaction Importance (CV k=5)",
              xlab = "",  
              horiz = TRUE,
              cex.names = 0.7,
              las = 1)

# Add error bars with adjusted bounds for horizontal orientation
segments(rev(lower_bound), bp, rev(upper_bound), bp, col = "red", lwd = 3)  # lwd = 3 for thicker lines

dev.off()
graphics.off() # reset margins

## Extract variables at the elbow for lm()
library(inflection)
rev(sort(var_importance_man$avg_var_props))

# Example variable importance data
bart_vars <- rev(sort(var_importance_man$avg_var_props))

# Step 1: Calculate the elbow point
set.seed(888)
elbow_index <- inflection::uik(x = seq_along(bart_vars), y = bart_vars)

# Step 2: Extract variables before the elbow
bart_vars <- names(bart_vars)[1:elbow_index]

# Mapping T1, T2, T3 to actual variables
variable_map <- list(
  T1 = "T_Newtech",
  T2 = "T_Specific",
  T3 = "T_Political"
)

# Replace T1, T2, T3 with their corresponding variables
bart_vars <- sapply(bart_vars, function(var) {
  for (key in names(variable_map)) {
    var <- sub(key, variable_map[[key]], var)
  }
  gsub(" x ", " * ", var)  # Replace ' x ' with ' * ' for interactions
})

## Check credible intervals for the marginal effects
bart_vars <- rev(sort(var_importance_man$avg_var_props))
bart_vars <- names(bart_vars)[names(bart_vars) %in% colnames(X_train)]


# Dynamically create interaction terms in X_train
for (var in names(bart_vars)) {
  if (!var %in% colnames(X_train)) {
    # Split the interaction term into components
    components <- strsplit(var, " \\* ")[[1]]
    if (all(components %in% colnames(X_train))) {
      # Compute the interaction term as the product of its components
      X_train[[var]] <- X_train[[components[1]]] * X_train[[components[2]]]
    } else {
      cat("Components not found for interaction term:", var, "\n")
    }
  }
}

## Plot effects
results <- lapply(bart_vars, function(var) {
  if (!var %in% colnames(X_train)) return(NULL)
  
  newdata_high <- X_train
  newdata_high[[var]] <- quantile(X_train[[var]], 0.9, na.rm = TRUE)
  newdata_low  <- X_train
  newdata_low[[var]]  <- quantile(X_train[[var]], 0.1, na.rm = TRUE)
  
  set.seed(888)
  pred_high <- bartMachine::bart_machine_get_posterior(bart_model_cv_man, newdata_high)$y_hat_posterior_samples
  set.seed(888)
  pred_low  <- bartMachine::bart_machine_get_posterior(bart_model_cv_man, newdata_low)$y_hat_posterior_samples
  
  posterior_diff <- pred_high - pred_low
  posterior_diff_mean <- colMeans(posterior_diff)
  
  data.frame(
    Variable = var,
    Lower = quantile(posterior_diff_mean, 0.025),
    Mean  = mean(posterior_diff_mean),
    Upper = quantile(posterior_diff_mean, 0.975)
  )
})

results_df <- do.call(rbind, results)

# filter based on a credible effect of zero - only one
results_df[results_df$Lower * results_df$Upper > 0, ]

# take top 10 for illustration
results_df$abs_mean <- abs(results_df$Mean)
results_df <- results_df[order(-results_df$abs_mean), ]
top10_df <- head(results_df, 10)


# Variables with credible evidence of an effect
print(top10_df)

xmax <- max(top10_df$Upper) + 0.025
xmin <- min(top10_df$Lower) - 0.025


################################################
### Figure A.9: Posterior Credible Intervals ###
################################################

png(filename="img/BART_PosteriorIntervals.png", width=1200, height=600, res=170)
ggplot(top10_df, aes(x = Variable, y = Mean, ymax=Upper, ymin=Lower)) + 
  geom_point(size = 3, alpha=.5) + 
  geom_errorbar(width = 0) + 
  coord_flip() +  # Flip coordinates for better readability 
  labs(title = "Posterior Credible Intervals for Interaction Variables (95% CI)",
       x = "", y = "") + 
  scale_y_continuous(
    breaks = seq(from = min(0, xmin), to = max(0, xmax), by = 0.5),
    limits = c(min(0, xmin), max(0, xmax))
  ) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 0.75) + 
  theme_minimal() 
dev.off()



########## DIAGNOSTICS ##########

## Extract the specific row with the minimum out-of-sample error
optimal_config <- bart_model_cv_man$cv_stats[which.min(bart_model_cv_man$cv_stats[, "oos_error"]), ]
m <- as.numeric(optimal_config["num_trees"])

# ....Or the lowest overall permutations
cv_stats_df <- as.data.frame(bart_model_cv_man$cv_stats)
summary_stats <- cv_stats_df %>%
  group_by(num_trees) %>%
  summarise(min_rmse = min(oos_error), mean_rmse = mean(oos_error), .groups = 'drop')

ggplot(cv_stats_df, aes(x = num_trees, y = oos_error, color = factor(nu), shape = factor(q))) +
  geom_point() +
  geom_line(aes(group = interaction(nu, q))) +
  labs(title = "RMSE by Number of Trees Across Configurations",
       x = "Number of Trees",
       y = "Out-of-Sample RMSE",
       color = "Degrees of Freedom",
       shape = "Quantile Level") +
  theme_minimal()

# Find the overall optimal number of trees not just by the lowest RMSE but also considering consistency
m <- as.numeric(summary_stats[which.min(summary_stats$min_rmse), "num_trees"])

# Additionally, explore which configurations yield the best performance at this optimal tree count
cv_stats_df %>%
  filter(num_trees == m) %>%
  arrange(oos_error)


## R^2 sensitivity
rsq1 <- k_fold_cv(as.data.frame(X_test), Y_test, k_folds = 10, use_missing_data=TRUE); rsq1

bart_model_cv_man$PseudoRsq
rsq1$PseudoRsq 

bart_model_cv_man$rmse 
rsq1$rmse


### RMSE via posterior prediction
train_preds_bart <- bart_model_cv_man$y_hat_train 
test_preds_bart <- predict(bart_model_cv_man, X_test) 
plot(Y_train, Y_train - train_preds_bart, main = "Residuals on Training Data")
plot(Y_test, Y_test - test_preds_bart, main = "Residuals on Test Data")
sqrt(mean((Y_train - train_preds_bart)^2)) # rmse for training set
sqrt(mean((Y_test - test_preds_bart)^2)) # rmse for test set






################################################
############# A.4.1: Topic Models ############## 
################################################

rm(list=setdiff(ls(), "df"))
set_bart_machine_num_cores(1) 
gc(verbose=TRUE, full=TRUE);gc(verbose=TRUE, full=TRUE);gc(verbose=TRUE, full=TRUE)
XLConnect::xlcFreeMemory()

# devtools::install_github("kbenoit/quanteda.dictionaries") 
library(quanteda)
library(quanteda.dictionaries)
library(stm)

## Simplifying spelling inconsistencies throughout the essays.
corp <- subset(df, !is.na(Y_Essay) & !is.na(Y_Likert))
corp <- corp %>% mutate(Text = stringi::stri_trans_general(Y_Essay, "Latin-ASCII")) # convert accented characters 
corp$Text <- gsub("#", "", corp$Text, ignore.case = T) # delete hashtags
corp$Text <- gsub("\\b\\d+(?!(,|º| degrees|km))\\b", "", corp$Text, perl = T) # clean up numbers unless unit of measurement
corp$Text <- gsub("xms|snx|nxs", "XNS", corp$Text, ignore.case = T) # fix obvious typo
corp$Text <- gsub("\\b(o\\.c\\.|oc|Orbital Convention)\\b", "Convention", corp$Text, ignore.case = T) # fix typo
corp$Text <- gsub("focuss", "focus", corp$Text, ignore.case = T) # fix obvious typo
corp$Text <- gsub("breif", "brief", corp$Text, ignore.case = T) # fix obvious typo
corp$Text <- gsub("acknolwedg", "acknowledg", corp$Text, ignore.case = T) # fix obvious typo
corp$Text <- gsub("centre", "center", corp$Text, ignore.case = T) # fix obvious typo
corp$Text <- gsub("baycent|barrycent|barclaycent|barcycent", "barycent", corp$Text, ignore.case = T) # fix obvious typo
corp$Text <- gsub("\\b(barycentr*)\\b", "barycenter", corp$Text, ignore.case = T) # fix obvious typo
corp$Text <- gsub("\\b(Article|Articles|Art|Arts|Art\\.|Arts\\.)\\b", "Article", corp$Text, ignore.case = T) # semantically identical
corp$Text <- gsub("articleii", "article ii", corp$Text, ignore.case = T) # fix obvious typo
corp$Text <- gsub("\\b(lissajous|lissajou|l)\\b", "Lissajous", corp$Text, ignore.case = T) # standardize references to L orbit
corp$Text <- gsub("center\\s+point", "centerpoint", corp$Text, ignore.case = T)
corp$Text <- tolower(corp$Text)


# Preprocessing
corp_full <- corpus(corp$Text, docvars=corp)



################################################
###### Figure A.10: Sensitivity Analysis #######
################################################

## devtools::install_github("matthewjdenny/preText")
library(preText)
preprocessed_documents <- factorial_preprocessing(corp_full, use_ngrams = TRUE, infrequent_term_threshold = 0.1, verbose = T)
preText_results <- preText(preprocessed_documents, dataset_name = "Corpus Sample", distance_method = "cosine", num_comparisons = 20, verbose = T)

head(preprocessed_documents$choices)
preText_score_plot(preText_results) # lower score => less risky

regression_coefficient_plot(preText_results, remove_intercept = TRUE)




## (1) Impact of T1 on emphasis | T2=1 specific --------
corp_spec <- subset(corp, T_Specific==1) # looking for evidence trying to protect principal

# Preprocess text
processed_spec <- textProcessor(corp_spec$Text, metadata = corp_spec,
                                stem = TRUE,
                                lowercase = TRUE,
                                removestopwords = TRUE,
                                removenumbers = TRUE,
                                removepunctuation = TRUE,
                                verbose = TRUE)

out_spec <- prepDocuments(processed_spec$documents, processed_spec$vocab, processed_spec$meta,
                          lower.thresh = 1, upper.thresh = Inf, verbose = TRUE) 


## (2) Impact of T2 on emphasis | T1=1 shock --------
corp_new <- subset(corp, T_Newtech==1) # looking for evidence trying to protect principal

# Preprocess text
processed_shock <- textProcessor(corp_new$Text, metadata = corp_new,
                                 stem = TRUE,
                                 lowercase = TRUE,
                                 removestopwords = TRUE,
                                 removenumbers = TRUE,
                                 removepunctuation = TRUE,
                                 verbose = TRUE)

out_shock <- prepDocuments(processed_shock$documents, processed_shock$vocab, processed_shock$meta,
                           lower.thresh = 1, upper.thresh = Inf, verbose = TRUE) 


## (3) Impact of T3 on emphasis --------
corp_for <- subset(corp, Y_Restrict>0) # looking for evidence trying to protect principal

# Preprocess text
processed_for <- textProcessor(corp_for$Text, metadata = corp_for,
                               stem = TRUE,
                               lowercase = TRUE,
                               removestopwords = TRUE,
                               removenumbers = TRUE,
                               removepunctuation = TRUE,
                               verbose = TRUE)

out_for <- prepDocuments(processed_for$documents, processed_for$vocab, processed_for$meta,
                         lower.thresh = 1, upper.thresh = Inf, verbose = TRUE)


################################################
####### Figure A.11: Topic Optimization ########
################################################

kResult_spec <- searchK(
  out_spec$documents, out_spec$vocab, data = out_spec$meta, 
  prevalence = ~ T_Newtech * T_Political, 
  K = c(6:20), 
  heldout.seed = 888, cores = 11, verbose = TRUE
)


graphics.off(); png(filename="stm_k_spec.png")
plot(kResult_spec)
dev.off()
k.spec <- 9


kResult_shock <- searchK(
  out_shock$documents, out_shock$vocab, data = out_shock$meta, 
  prevalence = ~ T_Specific * T_Political + X_Pre_JWT, 
  K = c(6:20), 
  heldout.seed = 888, cores = 11, verbose = TRUE
)


graphics.off(); png(filename="stm_k_shock.png")
plot(kResult_shock)
dev.off()
k.shock <- 9



kResult_for <- searchK(
  out_for$documents, out_for$vocab, data = out_for$meta, 
  prevalence = ~ T_Newtech * T_Specific * T_Political, 
  K = c(6:20), 
  heldout.seed = 888, cores = 11, verbose = TRUE
)

graphics.off(); png(filename="stm_k_for.png")
plot(kResult_for)
dev.off()
k.for <- 12



################################################
########## Figure A.12: Estimate STMs ##########
################################################

## Fit STM 1 with seed words
fit_spec <- stm(documents = out_spec$documents, vocab = out_spec$vocab, data = out_spec$meta, 
                prevalence = ~ T_Newtech + T_Political, 
                K = k.spec,
                init.type = "Spectral")

## Check topics
stm_spec <- estimateEffect(1:k.spec ~ T_Newtech + T_Political, 
                           stmobj = fit_spec,
                           meta = out_spec$meta,
                           uncertainty = "Global")

## Label topics using high-frequency words
labelTopics(fit_spec)
topic_labels_spec <- c("* Perceived Ordinary", # odd-numb, predict, satellit, criterion, creator
                       "Reasoning", # fail, arguendo, judgment, lissaj, met
                       "Case Facts", # ellip, curv, cycl, centerpoint, format
                       "Obligation (A4)", # forc, ratific, right, sovereign, withdraw
                       "* Clearly Restricted", # worldwid, lobbi, betray, enforc, countri 
                       "* Perceived Novelty", # lissajou, sun, moon, relat, miss
                       "Framer Intentions", # overcrowd, danger, preambl, damag, rational
                       "* Clearly Unrestricted", # lissajou, fail, indirect, baryc, experience
                       "Treaty Spirit" # discretion, plain meani, motiv, spirit, express 
)

## Plot STM 1
graphics.off()  # Clear all graphics devices
png("STM_FXPlots_T1.png", width = 1000, height = 1000, res = 230)
par(mfrow = c(1, 1))  # Set layout to 5 rows, 1 column

plot.estimateEffect(stm_spec, 
                    covariate = "T_Newtech", 
                    topics = 1:k.spec,
                    model = fit_spec,
                    method = "difference",
                    cov.value1 = 1, # treatment
                    cov.value2 = 0, # control
                    xlim = c(-1.3, .3),
                    ci.level=0.8, # (0.83 => p=0.05)
                    labeltype = "custom",
                    custom.labels = topic_labels_spec,
                    main = "Causes of Rule Disprution:\nInterpretations Under Specificity\n(T2=1 Fixed)",
                    xlab = "Mundane Tech (T1=0)  vs  Novel Tech (T1=1)"
)



## Fit STM 2 with seed words
fit_shock <- stm(documents = out_shock$documents, vocab = out_shock$vocab, data = out_shock$meta, 
                 prevalence = ~ T_Specific * T_Political + X_Pre_JWT,
                 K = k.shock,
                 init.type = "Spectral")

## Check topics
stm_shock <- estimateEffect(1:k.shock ~ T_Specific * T_Political + X_Pre_JWT,
                            stmobj = fit_shock,
                            meta = out_shock$meta,
                            uncertainty = "Global")


## Label topics using high-frequency words
labelTopics(fit_shock)
topic_labels_shock <- c("Policy Implications", #  debri, multipl, earth, collision, knot
                        "* Restricted by Analogy", # restrict, high, low, moder, undoubt, past
                        "* Law Permits", # foreclos, impos, prong, permit, criterion
                        "Sets Precedent", # advanc, futur, innov, floodgat, legisl, potential
                        "* Tech Excluded", #  fail, meet, less, facial, element
                        "* Update Needed", # interv, predict, profit, invest, update
                        "Obligation (A4)", # review, countri, law, withdraw, member, agreement
                        "* Fair Interpretation", # persuad, sound, scheme, obviat, certain
                        "* Framer Intentions" # match, deduc, region, fruit, framer, celesti, broad
)

## Plot STM 2
graphics.off()  # Clear all graphics devices
png("STM_FXPlots_T2.png", width = 1000, height = 1000, res = 230)
par(mfrow = c(1, 1))  # Set layout to 5 rows, 1 column

plot.estimateEffect(stm_shock, 
                    covariate = "T_Specific", 
                    topics = 1:k.shock,
                    model = fit_shock,
                    method = "difference",
                    cov.value1 = 0, # control
                    cov.value2 = 1, # treatment
                    xlim = c(-1.5, .6), 
                    ci.level=0.8, # (0.83 => p=0.05)
                    labeltype = "custom",
                    custom.labels = topic_labels_shock,
                    main = "Generality and Convergence:\nInterpretations About Novelty\n(T1=1 Fixed)",
                    xlab = "Specific (T2=1)     vs     General (T2=0)"
)




## Fit STM with seed words
fit_for <- stm(documents = out_for$documents, vocab = out_for$vocab, data = out_for$meta, 
               prevalence = ~ T_Newtech * T_Specific * T_Political, 
               K = k.for,
               init.type = "Spectral")

## Check topics
stm_for <- estimateEffect(1:k.for ~ T_Newtech * T_Specific * T_Political,
                          stmobj = fit_for,
                          meta = out_for$meta,
                          uncertainty = "Global")

## Label topics using high-frequency words
labelTopics(fit_for)
topic_labels_for <- c("Obligations", # purview, restrict, convent, articl, appli
                      "Framer Intentions", # framer, stifl, drafter, reinterpret, tradit
                      "* Reasoning", # high, low, median, travers, degre
                      "* Reputation", # worldwid, stark, countri, option, diplomat, secur
                      "Novel Tech", # lissajou, centerpoint, graviti, sun, moon
                      "* Credibility", # contest, logic, proof, wrinkl, defer 
                      "Legitimacy", # legitimacy, agreement, limit, enforc, opinion
                      "Analogies", # equival, belong, measur, third, satellit
                      "Treaty Spirit", # region, fruit, safeti, promot, enjoy
                      "Render Advice", # maintain, regular, revolut, convent, arti, sought 
                      "* Plain Meaning", # plain-meaning, disappoint, rigid, convent, prong
                      "Mundane Tech") # repetit, unchang, predict, flight, eccent

## Plot STM 3
graphics.off()  # Clear all graphics devices
png("STM_FXPlots_T3.png", width = 1000, height = 1000, res = 230)
par(mfrow = c(1, 1))  # Set layout to 5 rows, 1 column

plot.estimateEffect(stm_for, 
                    covariate = "T_Political", 
                    topics = 1:k.for,
                    model = fit_for,
                    method = "difference",
                    cov.value1 = 1, # treatment
                    cov.value2 = 0, # control
                    xlim = c(-.6, .2),
                    ci.level=0.8, # (0.83 => p=0.05)
                    moderator = "T_Specific",
                    moderator.value = 0,
                    labeltype = "custom",
                    custom.labels = topic_labels_for,
                    main = "Reasons to Restrict:\nPersuading the Principal\n (Y > 0 Fixed)",
                    xlab = "Impartial (T3=0)     vs     Politicized (T3=1)"
)
dev.off()


## Model validation (weak correlations imply high distinctiveness)
plot(topicCorr(fit_for))
plot(topicCorr(fit_for, method="huge", verbose = TRUE))
plot(topicCorr(fit_shock))
plot(topicCorr(fit_shock, method="huge", verbose = TRUE))
plot(topicCorr(fit_spec))
plot(topicCorr(fit_spec, method="huge", verbose = TRUE))


######### Topic proportions ######### 

## Extract theta from the models
word_df <- data.frame(topic_labels_shock)
proportion <- as.data.frame(colSums(fit_shock$theta/nrow(fit_shock$theta)))
word_df <- cbind(word_df, proportion)
colnames(word_df) <- c("Labels", "Probability")

## Sort the dataframe
word_df <- word_df %>% arrange(-Probability)
word_df$Labels <- factor(word_df$Labels, levels = rev(word_df$Labels))
word_df$Probability <- as.numeric(word_df$Probability)
word_df$Probability <- round(word_df$Probability, 4)

## Plot 
ggplot(word_df, aes(x = Labels, y = Probability)) + 
  geom_bar(stat = "identity") + 
  scale_y_continuous(breaks = c(0, 0.25), limits = c(0, 0.25), expand = c(0, 0)) + #change breaks and limits as you need
  coord_flip() + 
  geom_text(aes(label = scales::percent(Probability)), #Scale in percent
            hjust = -0.25, size = 4,
            position = position_dodge(width = 1),
            inherit.aes = TRUE) + 
  theme(panel.border = element_blank()) +
  theme_minimal() +
  labs(title = "STM Model 1 (Impact of Tech Change on T2=1)",
       x = "", y = "") 


## Extract theta from the stm-model
word_df <- data.frame(topic_labels_spec)
proportion <- as.data.frame(colSums(fit_spec$theta/nrow(fit_spec$theta)))
word_df <- cbind(word_df, proportion)
colnames(word_df) <- c("Labels", "Probability")

## Sort the dataframe
word_df <- word_df %>% arrange(-Probability)
word_df$Labels <- factor(word_df$Labels, levels = rev(word_df$Labels))
word_df$Probability <- as.numeric(word_df$Probability)
word_df$Probability <- round(word_df$Probability, 4)

## Plot graph
ggplot(word_df, aes(x = Labels, y = Probability)) + 
  geom_bar(stat = "identity") + 
  scale_y_continuous(breaks = c(0, 0.25), limits = c(0, 0.25), expand = c(0, 0)) + #change breaks and limits as you need
  coord_flip() + 
  geom_text(aes(label = scales::percent(Probability)), #Scale in percent
            hjust = -0.25, size = 4,
            position = position_dodge(width = 1),
            inherit.aes = TRUE) + 
  theme(panel.border = element_blank()) +
  theme_minimal() +
  labs(title = "STM Model 2 (Rule Resilience Given T1=1)",
       x = "", y = "") 


## Extract theta from the stm-model
word_df <- data.frame(topic_labels_for)
proportion <- as.data.frame(colSums(fit_for$theta/nrow(fit_for$theta)))
word_df <- cbind(word_df, proportion)
colnames(word_df) <- c("Labels", "Probability")

## Sort the dataframe
word_df <- word_df %>% arrange(-Probability)
word_df$Labels <- factor(word_df$Labels, levels = rev(word_df$Labels))
word_df$Probability <- as.numeric(word_df$Probability)
word_df$Probability <- round(word_df$Probability, 4)

## Plot graph
ggplot(word_df, aes(x = Labels, y = Probability)) + 
  geom_bar(stat = "identity") + 
  scale_y_continuous(breaks = c(0, 0.25), limits = c(0, 0.25), expand = c(0, 0)) + #change breaks and limits as you need
  coord_flip() + 
  geom_text(aes(label = scales::percent(Probability)), #Scale in percent
            hjust = -0.25, size = 4,
            position = position_dodge(width = 1),
            inherit.aes = TRUE) + 
  theme(panel.border = element_blank()) +
  theme_minimal() +
  labs(title = "STM Model 3 (Client Advice for Y > 0)",
       x = "", y = "") 


################################################
######### Table A.5: Response Excerpts ######### 
################################################


## KWIC Table
library(quanteda)

# Ensure the essays are treated as individual documents
corpus_essays <- corpus(df$Y_Essay, docnames = as.character(1:nrow(df)))

# Columns that indicate treatment conditions
treatment_cols <- grep("^T\\d{3}$", names(df), value = TRUE)

# Create a vector to store the treatment condition for each row
df$Z_Condition <- apply(df[, treatment_cols], 1, function(x) names(x)[x == 1])

# Generate KWIC output
kwic_output <- kwic(tokens(corpus_essays), pattern = "credib*|reputation|alleg*|noncompl*|comply|complia*|overreach|overrule|bias|justif", window = 20)

# Convert kwic_output to DataFrame
kwic_df <- data.frame(
  docname = kwic_output$docname,
  pre = kwic_output$pre,
  keyword = kwic_output$keyword,
  post = kwic_output$post,
  stringsAsFactors = FALSE
)

# Concatenate pre, keyword, and post
kwic_df$`KWIC Window` <- paste(kwic_df$pre, kwic_df$keyword, kwic_df$post)

# Append Synthetic_ID using the docname to match
kwic_df$`Participant ID` <- df$Synthetic_ID[as.numeric(kwic_df$docname)]
kwic_df$`Z Condition` <- df$Z_Condition[as.numeric(kwic_df$docname)]

# Select relevant cols only
kwic_df <- kwic_df[, c("KWIC Window", "Participant ID", "Z Condition")]

# Convert DataFrame to LaTeX table
library(xtable)
latex_table <- xtable(kwic_df)

# Create LaTeX table with tabularx
print(latex_table, include.rownames = FALSE, floating = FALSE, hline.after = NULL, 
      sanitize.text.function = function(x){x}, 
      tabular.environment = "tabularx", 
      width = "\\textwidth", 
      align = "p{7cm}cc") # Adjust the alignment as needed, 'X' for automatic line breaks



################################################
####### Figure A.13: Mediation Analysis ######## 
################################################

library(mediation)
library(dplyr)
library(tidyr)


treatments <- c("T_Newtech", "T_Specific", "T_Political")
mediators  <- c("Y_CiteAmbig", "Y_CiteArticles1", "Y_CiteArticles2", "Y_CiteArticles3",
                "Y_CiteArticles4", "Y_CitePrecedent", "Y_CitePreamble", "Y_CiteNovelty",
                "Y_CiteEffects", "Y_CiteDefinitions", "Y_CiteFears", "Y_Nchar",
                "Timer_Read", "Timer_Write", "Timer_Clicks")
outcome    <- "Y_Likert"

results <- list()

library(mediation)

mediation_results <- list()

for (med in mediators) {
  message("\n=== Mediator: ", med, " ===")
  res_list <- list()
  
  # Mediator and outcome models include all treatment terms
  m.form <- as.formula(paste(med, "~ T_Newtech + T_Specific + T_Political + T_Newtech * T_Specific"))
  y.form <- as.formula(paste(outcome, "~ T_Newtech + T_Specific + T_Political + T_Newtech * T_Specific +", med))
  
  m.mod <- lm(m.form, data = df)
  y.mod <- lm(y.form, data = df)
  
  for (treat in treatments) {
    message("  →  ", treat)
    res <- tryCatch(
      mediate(
        model.m = m.mod,
        model.y = y.mod,
        treat = treat,
        mediator = med,
        boot = TRUE,
        sims = 1000
      ),
      error = function(e) { message("     ❌ ", e$message); NULL }
    )
    if (!is.null(res)) res_list[[treat]] <- res
  }
  mediation_results[[med]] <- res_list
}

# extract numeric summaries
med_table <- bind_rows(lapply(names(results), function(med) {
  bind_rows(lapply(names(results[[med]]), function(treat) {
    s <- summary(results[[med]][[treat]])
    data.frame(
      Mediator = med,
      Treatment = treat,
      ACME = s$d0, ACME_L = s$d0.ci[1], ACME_U = s$d0.ci[2],
      ADE  = s$z0, ADE_L  = s$z0.ci[1], ADE_U  = s$z0.ci[2],
      Total = s$tau.coef, Total_L = s$tau.ci[1], Total_U = s$tau.ci[2],
      PropMed = s$n0, Prop_L = s$n0.ci[1], Prop_U = s$n0.ci[2]
    )
  }))
}))



# prepare empty data frame
library(dplyr)

# ---- 1️⃣ Collect all mediation results into a table ----
mediation_table <- data.frame()

for (med in names(mediation_results)) {
  for (treat in names(mediation_results[[med]])) {
    med_obj <- mediation_results[[med]][[treat]]
    if (!is.null(med_obj)) {
      s <- summary(med_obj)
      mediation_table <- bind_rows(
        mediation_table,
        data.frame(
          Mediator      = med,
          Treatment     = treat,
          ACME          = s$d0,
          ACME_LowerCI  = s$d0.ci[1],
          ACME_UpperCI  = s$d0.ci[2],
          ADE           = s$z0,
          ADE_LowerCI   = s$z0.ci[1],
          ADE_UpperCI   = s$z0.ci[2],
          TotalEffect   = s$tau.coef,
          Total_LowerCI = s$tau.ci[1],
          Total_UpperCI = s$tau.ci[2],
          PropMediated  = s$n0,
          Prop_LowerCI  = s$n0.ci[1],
          Prop_UpperCI  = s$n0.ci[2]
        )
      )
    }
  }
}

# --- filter grid to only mediator × treatment pairs that exist ---
# --- identify significant ACME effects ---------------------------------------

sig_pairs <- mediation_table %>%
  filter(!is.na(ACME_LowerCI), !is.na(ACME_UpperCI)) %>%
  filter(ACME_LowerCI * ACME_UpperCI > 0) %>%
  distinct(Mediator, Treatment)

mediators_sig <- unique(sig_pairs$Mediator)


# safety check
if (length(mediators_sig) == 0) {
  message("No mediators with significant ACME — nothing to plot.")
  quit(save = "no")  # or simply return() if inside a function
}


meds <- unique(sig_pairs$Mediator)
treats <- unique(sig_pairs$Treatment)

nrow <- length(meds)
ncol <- length(treats)

px_per_panel <- 900
dpi <- 300

png("Mediation_sigACME_grid.png",
    width  = ncol * px_per_panel,
    height = nrow * px_per_panel,
    res    = dpi)

par(mfrow = c(nrow, ncol),
    mar = c(2, 4, 3, 1),
    oma = c(3, 4, 2, 1),
    cex = 0.9)

for (med in meds) {
  for (treat in treats) {
    
    # skip if combination not in results
    if (!med %in% names(mediation_results)) { plot.new(); next }
    if (!treat %in% names(mediation_results[[med]])) { plot.new(); next }
    
    med_obj <- mediation_results[[med]][[treat]]
    if (is.null(med_obj)) { plot.new(); next }
    
    s <- summary(med_obj)
    xlo <- min(s$d0.ci[1], s$z0.ci[1], s$tau.ci[1])
    xhi <- max(s$d0.ci[2], s$z0.ci[2], s$tau.ci[2])
    pad <- max(0.5, 0.15 * (xhi - xlo))
    xlim_use <- c(xlo - pad, xhi + pad)
    
    plot(
      med_obj,
      main = paste0("Mediator: ", med, "\nTreatment: ", treat),
      xlim = xlim_use,
      cex.main = 0.9
    )
  }
}

# add outer labels
#mtext("Treatments", side = 3, outer = TRUE, line = 0.5, cex = 1)
#mtext("Mediators", side = 2, outer = TRUE, line = 2, cex = 1)

dev.off()
message("✅ Saved Mediation_sigACME_grid.png with ", nrow, " mediators × ", ncol, " treatments.")



# Focused plot
sig_df <- mediation_table %>%
  filter(!is.na(ACME_LowerCI), !is.na(ACME_UpperCI)) %>%
  mutate(Significant = ACME_LowerCI * ACME_UpperCI > 0)

# ACME estimates with 95% CIs by mediator and treatment
png(filename="Mediation_signifACME.png", width=1200, height=600, res=150)
ggplot(sig_df, aes(x = Treatment, y = ACME, color = Significant)) +
  geom_pointrange(aes(ymin = ACME_LowerCI, ymax = ACME_UpperCI)) +
  facet_wrap(~ Mediator, scales = "free_y") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(y = "Average Causal Mediation Effect (ACME)",
       x = "",
       title = "")
dev.off()




